mydata <- read.csv2("Apartments.csv") 
head(mydata)
##   ID PriceM2 Size Distance Floor Parking Age
## 1  1    2.64   70      8.4     3       0   3
## 2  2    3.80   51      0.6     1       1   3
## 3  3    2.66   68      8.4     3       0   3
## 4  4    2.85   70      8.1     3       0   3
## 5  5    2.64   68      5.4     3       1   3
## 6  6    2.77   76      3.6     2       0   3

Description of the variables:

mydata$ParkingF <- factor(mydata$Parking, 
                          levels = c(0, 1),
                          labels = c("No", "Yes"))

mydata$AgeF <- factor(mydata$Age, 
                            levels = c(1, 2, 3),
                            labels = c("New", "≤10 yrs", ">10 yrs"))

summary(mydata[c(-1, -6, -7)])
##     PriceM2           Size          Distance         Floor      
##  Min.   :2.400   Min.   :40.00   Min.   :0.100   Min.   :0.000  
##  1st Qu.:2.808   1st Qu.:57.00   1st Qu.:1.125   1st Qu.:2.000  
##  Median :3.250   Median :63.00   Median :3.600   Median :3.000  
##  Mean   :3.434   Mean   :63.01   Mean   :3.761   Mean   :3.634  
##  3rd Qu.:3.915   3rd Qu.:70.00   3rd Qu.:5.575   3rd Qu.:5.000  
##  Max.   :4.980   Max.   :82.00   Max.   :8.700   Max.   :9.000  
##  ParkingF      AgeF   
##  No :29   New    :48  
##  Yes:53   ≤10 yrs:20  
##           >10 yrs:14  
##                       
##                       
## 
library(GGally)  # Load tools for scatterplot matrices

ggpairs(
  mydata[c("PriceM2", "Size", "Distance", "Floor", "ParkingF", "AgeF")],  # Variables to display
  upper = list(
    continuous = wrap("smooth", method = "lm", se = FALSE, alpha = 0.6)  # Scatterplots with OLS line
  ),
  lower = list(continuous = "cor"),   # Pearson correlation coefficients
  diag  = list(continuous = "densityDiag")  # Density plots on the diagonal
) +
theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 7),
    axis.text.y = element_text(size = 7)
)

fit <- lm(PriceM2 ~ Size + Distance + Floor + ParkingF + AgeF, 
           data = mydata)
library(car)
## Loading required package: carData
vif(fit) #Checking multicolinearity
##              GVIF Df GVIF^(1/(2*Df))
## Size     1.788840  1        1.337475
## Distance 1.701312  1        1.304344
## Floor    1.200754  1        1.095789
## ParkingF 1.313215  1        1.145956
## AgeF     1.241904  2        1.055655
mean(vif(fit)[, "GVIF^(1/(2*Df))"])
## [1] 1.187844
mydata$StdResid <- rstandard(fit)
hist(mydata$StdResid, 
     xlab = "Standardised residuals", 
     ylab = "Frequency", 
     main = "Histogram of standardised residuals")

shapiro.test(mydata$StdResid)
## 
##  Shapiro-Wilk normality test
## 
## data:  mydata$StdResid
## W = 0.98681, p-value = 0.5686
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Calculate Cook's distance
mydata$CooksD <- cooks.distance(fit)

# Define Cook's distance threshold
cd_threshold <- 4/nrow(mydata)

# Display apartments with Cook's distance above the threshold
influential_units <- mydata %>%
  filter(CooksD > cd_threshold) %>%
  arrange(desc(CooksD)) %>%
  select(ID, PriceM2, Size, Distance, Floor, ParkingF, AgeF, CooksD)

influential_units %>%
  mutate(CooksD = round(CooksD, 3))
##   ID PriceM2 Size Distance Floor ParkingF    AgeF CooksD
## 1 81    4.90   66      0.1     4      Yes     New  0.075
## 2 75    4.60   68      1.2     4      Yes     New  0.051
## 3 10    3.57   59      1.2     8      Yes >10 yrs  0.051
## 4  4    2.85   70      8.1     3       No >10 yrs  0.049
mydata$StdFittedValues <- scale(fit$fitted.values)

library(car)
scatterplot(y = mydata$StdResid, x = mydata$StdFittedValues,
            ylab = "Standardised residuals",
            xlab = "Standardised fitted values",
            boxplots = FALSE,
            regLine = FALSE,
            smooth = FALSE)

library(olsrr)
## 
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
## 
##     rivers
ols_test_breusch_pagan(fit)
## 
##  Breusch Pagan Test for Heteroskedasticity
##  -----------------------------------------
##  Ho: the variance is constant            
##  Ha: the variance is not constant        
## 
##                Data                 
##  -----------------------------------
##  Response : PriceM2 
##  Variables: fitted values of PriceM2 
## 
##         Test Summary         
##  ----------------------------
##  DF            =    1 
##  Chi2          =    2.069615 
##  Prob > Chi2   =    0.1502595
summary(fit) 
## 
## Call:
## lm(formula = PriceM2 ~ Size + Distance + Floor + ParkingF + AgeF, 
##     data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73396 -0.22692 -0.02825  0.22486  0.87658 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.16921    0.36885  14.015  < 2e-16 ***
## Size        -0.02511    0.00556  -4.517 2.29e-05 ***
## Distance    -0.12199    0.01835  -6.648 4.21e-09 ***
## Floor        0.04792    0.02197   2.181 0.032285 *  
## ParkingFYes  0.33227    0.08990   3.696 0.000414 ***
## AgeF≤10 yrs -0.14155    0.09720  -1.456 0.149484    
## AgeF>10 yrs -0.28479    0.10999  -2.589 0.011551 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3396 on 75 degrees of freedom
## Multiple R-squared:  0.807,  Adjusted R-squared:  0.7916 
## F-statistic: 52.27 on 6 and 75 DF,  p-value: < 2.2e-16
sqrt(summary(fit)$r.squared) 
## [1] 0.8983382
library(lm.beta)
lm.beta(fit)
## 
## Call:
## lm(formula = PriceM2 ~ Size + Distance + Floor + ParkingF + AgeF, 
##     data = mydata)
## 
## Standardized Coefficients::
## (Intercept)        Size    Distance       Floor ParkingFYes 
##          NA -0.30645540 -0.43985221  0.12125659  0.21485378 
## AgeF≤10 yrs AgeF>10 yrs 
## -0.08221483 -0.14492675
fitNEW <- lm(
  PriceM2 ~ Size + I(Size^2) +
             Distance + I(Distance^2) +
             Floor + ParkingF + AgeF,
  data = mydata
)

mydata$StdResid <- round(rstandard(fitNEW), 3) 

mydata$StdFittedValues <- scale(fitNEW$fitted.values)

library(car)
scatterplot(y = mydata$StdResid, x = mydata$StdFittedValues,
            ylab = "Standardized residuals",
            xlab = "Standardized fitted values",
            boxplots = FALSE,
            regLine = FALSE,
            smooth = FALSE)

summary(fitNEW)
## 
## Call:
## lm(formula = PriceM2 ~ Size + I(Size^2) + Distance + I(Distance^2) + 
##     Floor + ParkingF + AgeF, data = mydata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8072 -0.1295  0.0215  0.1549  0.6680 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.6265496  1.2696166   4.432 3.23e-05 ***
## Size          -0.0354422  0.0403222  -0.879  0.38230    
## I(Size^2)      0.0001386  0.0003156   0.439  0.66190    
## Distance      -0.3914330  0.0504730  -7.755 4.06e-11 ***
## I(Distance^2)  0.0305664  0.0054130   5.647 2.95e-07 ***
## Floor          0.0415768  0.0185466   2.242  0.02802 *  
## ParkingFYes    0.3342930  0.0755459   4.425 3.31e-05 ***
## AgeF≤10 yrs   -0.0472774  0.0833000  -0.568  0.57208    
## AgeF>10 yrs   -0.3167237  0.0928486  -3.411  0.00106 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2852 on 73 degrees of freedom
## Multiple R-squared:  0.8675, Adjusted R-squared:  0.853 
## F-statistic: 59.77 on 8 and 73 DF,  p-value: < 2.2e-16
mydata$FittedValues <- fitted.values(fit) 
mydata$Residuals <- residuals(fit) 
mydata[2 , colnames(mydata) %in% c("ID", "PriceM2", "FittedValues", "Residuals")]
##   ID PriceM2 FittedValues  Residuals
## 2  2     3.8     3.910591 -0.1105908